74 


ELSEVIER 


journal  homepage:  www.elsevier.com/locate/epilepsyres 


Epileptic  seizures  from  abnormal  networks:  Why 
some  seizures  defy  predictability 

William  S.  Anderson3  *,  Feraz  Azharb1,  Pawel  Kudelac’* 1 2, 

Gregory  K.  Bergeyc3,  Piotr  J.  Franaszczukd  e’4 

a  The  Johns  Hopkins  University  School  of  Medicine,  Department  of  Neurosurgery,  Meyer  5-109E,  600  North  Wolfe  Street, 
Baltimore,  MD  21287,  USA 

b  The  Johns  Hopkins  University  School  of  Medicine,  Department  of  Neurosurgery,  Meyer  5-157,  600  North  Wolfe  Street, 
Baltimore,  MD  21287,  USA 

c  The  Johns  Hopkins  University  School  of  Medicine,  Department  of  Neurology,  Meyer  2- 147,  600  North  Wolfe  Street,  Baltimore, 
MD  21287,  USA 

d  The  Johns  Hopkins  University  School  of  Medicine,  Department  of  Neurology,  USA 

e  U.S.  Army  Research  Laboratory,  Human  Research  and  Engineering  Directorate,  Aberdeen  Proving  Ground,  2800  Powder  Mill 
Road,  Adelphi,  MD  20783,  USA 

Received  15  June  2011;  received  in  revised  form  19  October  2011;  accepted  18  November  2011 
Available  online  12  December  2011 


KEYWORDS 

Computational 

simulation; 

Neural  network 
model; 

Seizure  prediction; 
Seizure  generation 


Summary  Seizure  prediction  has  proven  to  be  difficult  in  clinically  realistic  environments.  Is  it 
possible  that  fluctuations  in  cortical  firing  could  influence  the  onset  of  seizures  in  an  ictal  zone? 
To  test  this,  we  have  now  used  neural  network  simulations  in  a  computational  model  of  cortex 
having  a  total  of  65,536  neurons  with  intercellular  wiring  patterned  after  histological  data.  A 
spatially  distributed  Poisson  driven  background  input  representing  the  activity  of  neighboring 
cortex  affected  1%  of  the  neurons.  Gamma  distributions  were  fit  to  the  interbursting  phase 
intervals,  a  non-parametric  test  for  randomness  was  applied,  and  a  dynamical  systems  analysis 
was  performed  to  search  for  period-1  orbits  in  the  intervals.  The  non-parametric  analysis  sug¬ 
gests  that  intervals  are  being  drawn  at  random  from  their  underlying  joint  distribution  and  the 
dynamical  systems  analysis  is  consistent  with  a  nondeterministic  dynamical  interpretation  of 
the  generation  of  bursting  phases.  These  results  imply  that  in  a  region  of  cortex  with  abnormal 
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connectivity  analogous  to  a  seizure  focus,  it  is  possible  to  initiate  seizure  activity  with  fluctu¬ 
ations  of  input  from  the  surrounding  cortical  regions.  These  findings  suggest  one  possibility  for 
ictal  generation  from  abnormal  focal  epileptic  networks.  This  mechanism  additionally  could  help 
explain  the  difficulty  in  predicting  partial  seizures  in  some  patients. 

©  2011  Elsevier  B.V.  All  rights  reserved. 


Introduction 

Epileptic  seizures  are  brief,  episodic  phenomena.  Partial 
seizures,  the  most  common  seizure  type,  arise  from  focal 
brain  regions  (e.g.  temporal,  parietal)  (Niedermeyer,  2005). 

While  in  some  instances  there  may  be  an  identifiable 
cause  for  the  seizures  (e.g.  tumor,  cavernoma,  hippocam¬ 
pal  sclerosis),  in  other  instances  no  clear  pathology  is 
determined.  The  hallmark  of  an  epileptic  seizure  is  the 
involvement  of  local  or  regional  neural  networks;  repeti¬ 
tive  firing  of  a  single  neuron  does  not  produce  symptoms 
without  this  network  involvement.  What  causes  the  interic- 
tal  to  ictal  transition?  A  typical  partial  seizure  lasts  less  than 
2  min  plus  any  postictal  state  (Afra  et  al.,  2008).  Therefore, 
even  if  a  patient  has  very  frequent  seizures,  the  majority 
of  time  is  spent  in  the  interictal  state.  While  some  seizures 
can  be  provoked  or  are  more  likely  to  occur  under  certain 
situations  (e.g.  sleep  deprivation,  photic  stimulation),  the 
majority  of  seizures  appear  to  occur  spontaneously  without 
known  association  with  definable  influences. 

There  has  been  considerable  interest  in  seizure  predic¬ 
tion  in  recent  years.  Obviously  if  seizures  could  be  reliably 
predicted,  then  the  option  for  targeted  therapy  exists  (e.g. 
stimulation),  or  at  least  the  patient  could  remove  them¬ 
selves  from  potentially  dangerous  situations.  The  underlying 
hypothesis  for  seizure  prediction  is  that  there  are  changes 
in  cerebral  dynamics  that  may  precede  the  clinical  seizure 
by  minutes  to  hours  (reviewed  in  Sackellares,  2008).  These 
changes  may  be  local  (i.e.  near  the  seizure  focus)  or  remote. 
These  changes  are  not  apparent  with  visual  analysis  of  the 
EEG,  even  with  intracranial  recording  arrays.  Some  groups 
have  identified  high  frequency  activity  that  may  signal  the 
onset  of  neocortical  partial  seizures,  but  this  is  an  exam¬ 
ple  of  improved  seizure  detection,  not  prediction  (Worrell 
etal.,  2004,  2008;  Bragin  etal.,  2010).  Reliable  seizure  pred¬ 
ication  has  been  challenging  and  even  the  most  enthusiastic 
proponents  of  the  prediction  hypothesis  acknowledge  the 
difficulties  with  current  algorithms  (Lehnertz  et  al.,  2007; 
Mormann  et  al.,  2007;  Andrzejak  et  al.,  2009). 

Seizure  prediction  may  be  difficult  due  to  rapid  bistable 
state  changes  at  the  time  of  ictal  onset  in  the  neocortex 
(Suffczynski  et  al.,  2006;  Lopes  da  Silva  et  al.,  2003). 

The  mechanisms  underlying  a  bistable  state  change 
may  be  quite  different  between  primary  generalized  (e.g. 
absence)  and  partial  epileptic  seizures.  A  bistable  state 
change  may  be  more  applicable  to  these  primary  general¬ 
ized  seizures  which  have  abrupt  bilateral  cerebral  onset. 
In  this  paper,  a  different  possible  mechanism  is  presented 
under  which  seizure  prediction  would  be  difficult  in  some 
patients  with  focal  seizure  onset. 

Knowing,  as  we  do,  that  partial  seizures  are  a  reflec¬ 
tion  of  transient  abnormal  regional  network  activity,  it 
is  reasonable  to  postulate  that  these  seizures  in  at  least 
some  (perhaps  many)  patients  result  from  abnormal  neural 


networks  (e.g.  the  epileptogenic  zone)  (Jacobs  etal.,  2000). 
We  describe  here  a  model  of  the  epileptogenic  zone  where 
the  epileptic  focus  is  represented  by  an  abnormal  neu¬ 
ral  network  that  has  very  slightly  altered  connectivity  so 
that,  while  seizures  only  occur  infrequently,  they  can  be 
triggered  by  normal  background  activity  originating  from 
outside  the  epileptogenic  zone.  This  background  activity 
could  be  influenced  by  various  physiologic  factors  (e.g. 
sleep),  but  nevertheless  this  background  activity  would  not 
result  in  seizure  activity  in  the  non-epileptic  brain.  This 
does  not  discount  the  possibility  that  some  changes  in  neural 
network  synchrony  may  occur  in  the  "normal”  brain  since 
the  cumulative  lifetime  incidence  of  unprovoked  seizures 
approaches  4%  (Hauser  et  al.,  1993).  Often  these  seizures 
are  provoked  (e.g.  medications  and  alcohol)  and  less  than 
half  of  these  patients  have  recurrent  seizures.  The  life¬ 
time  cumulative  risk  of  developing  epilepsy  only  ranges 
from  1 .4%  to  3.3%  (Krumholz  et  al.,  2007;  Berg  and  Shinnar, 
1991).  In  this  model,  however,  where  normal  background 
activity,  occasionally  or  rarely  produces  a  seizure  in  abnor¬ 
mal  regional  networks,  seizure  prediction  would  be  difficult 
since  detectable  preictal  changes  would  not  be  present;  the 
first  changes  would  in  fact  be  seizure  initiation. 

Epileptic  networks  in  neocortex  or  the  hippocampus 
show  anatomical  changes  compared  to  normal  tissue  (Jacobs 
et  al.,  2000;  Sallin  et  al.,  1995).  These  changes  can  progress 
with  time  (Sallin  et  al.,  1995;  Arellano  et  al.,  2004).  This 
could  result  in  neuronal  networks  more  amenable  to  seizure 
generation  (electrical  or  clinical)  over  large  regional  areas. 
There  is  a  complex  interrelationship,  much  of  it  not  well 
understood,  between  neurons  which  are  dysfunctional  and 
the  neural  networks  which  can  promote  seizures  (Leussis 
and  Heinrichs,  2007;  Kumar  et  al.,  2007;  Swann  et  al., 
2007).  Even  in  the  non-epileptic  brain,  excitatory  connec¬ 
tions  predominate  with  80—90%  of  synapses  being  excitatory 
(Braitenberg  and  Schiiz,  1998). 

With  neuronal  network  simulations  it  is  possible  to  con¬ 
trol,  study,  and  quickly  change  the  various  influences  on 
network  behavior.  Recently,  we  presented  the  results  of 
computational  simulation  studies  examining  the  role  of 
external  field  stimulation  on  ongoing  bursting  activity  in  a 
neural  network  (Anderson  et  al.,  2007,  2009).  The  corti¬ 
cal  model  used  in  these  studies  consists  of  discrete  single 
compartment  Hodgkin— Huxley  type  cells  which  are  spa¬ 
tially  arranged  in  a  realistic  fashion  having  both  a  layered 
and  columnar  structure.  Since  neural  network  behavior 
reflects  the  aggregate  output  of  the  component  neurons, 
single  compartment  neurons  allow  greater  computational 
efficiency  and  the  ability  to  model  larger  networks  in  studies 
of  network  behavior.  Arrangements  of  connected  simu¬ 
lated  neurons  in  this  manner  can  demonstrate  spontaneous 
bursting  phases  and  have  spatial  characteristics  similar  to 
seizures  recorded  from  humans  (Anderson  et  al. ,  2007,  2009; 
Kudela  et  al.,  1997,  2003a, b,  2005;  Franaszczuk  et  al., 
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2003).  We  now  present  the  results  of  a  similar  neuronal 
network  model  with  random  surrounding  background  inputs. 
The  goal  of  this  study  was  to  investigate  the  statistical  struc¬ 
ture  of  the  resulting  bursting  network  activity  to  seek  the 
presence  or  absence  of  predictable  patterns  in  the  behavior. 

Materials  and  methods 
Computational  model  format 


background  inputs  were  the  same  used  in  the  cell  to  cell  connec¬ 
tions,  and  followed  the  same  rise  and  decay  times  appropriate  for 
postsynaptic  potentials. 

The  pseudo-random  number  generator  used  for  the  application 
of  the  noise  pulses  was  a  linear  congruential  generator  implemented 
with  the  C-function  drand48,  with  an  intrinsic  period  of  281  x  1012. 
For  a  30  s  simulation  and  1 0  5  s  time-step,  this  function  was  called 
1.966  x  109  times  for  1%  of  the  cells  undergoing  background  input. 
The  period  length  for  the  pseudo-random  generator  is  143,000  times 
larger  than  this  number. 


The  individual  neurons  in  this  neocortical  model  were  represented 
by  single  compartment  neurons  bearing  synaptic  connections  from 
the  rest  of  the  network,  and  embedded  with  a  fixed  set  of  ionic 
conductances.  The  membrane  potential  varies  as: 

dV 

Cm  ^  =  I  syn  —  ill  a  —  fca  —  /k  —  iR(Ca)  “  /a  —  h 

The  individual  currents  include  the  input  synaptic  current  /syn, 
inward  sodium  and  calcium  currents  /Na  and  /Ca,  outward  potas¬ 
sium  currents  including  a  delayed  rectifier  current  /«,  a  calcium 
dependent  potassium  current  /K(ca),  a  transient  potassium  current 
/A,  and  a  leakage  current  lL  (Av-Ron,  1994).  The  minicolumns  used 
in  the  simulation  consist  of  16  cells  with  intrinsic  intracolumnar 
wiring,  adapted  from  the  neocortical  work  of  Douglas  and  Martin 
(2004),  as  described  more  fully  in  previous  studies  (Anderson  et  al., 
2007,  2009).  This  is  both  for  its  ease  of  implementation  compu¬ 
tationally  and  for  its  experimental  support  in  somatosensory  and 
visual  cortex  (Douglas  and  Martin,  2004).  The  geometry  imposed  on 
a  computational  model  becomes  relevant  when  studying  any  spa¬ 
tially  dependent  effects  on  the  resultant  spreading  activity.  The 
minicolumns  in  this  simulation  have  a  25  pm  center-to-center  spac¬ 
ing  in  a  square  lattice  repeating  structure.  The  total  number  of 
cells  examined  was  65,536,  representing  a  simulated  cortical  sur¬ 
face  area  of  1.6  mm  x  1.6  mm.  Fig.  1  demonstrates  a  schematic  of 
the  intracolumnar  excitatory  cell  connections  and  the  organization 
of  the  minicolumns  in  planar  space  as  well  as  snapshots  of  the 
resultant  activity  in  the  layer  ll/lll  pyramidal  cell  component  dur¬ 
ing  model  bursting  activity.  The  model  connectivity  and  synaptic 
currents  are  described  further  in  the  Supplementary  material. 

The  base  connection  pattern  studied  in  this  report  is  represen¬ 
tative  of  one  that  can  produce  robust  bursting  as  previously  studied 
(Anderson  et  al.,  2007).  The  numbers  of  extra-columnar  connec¬ 
tions  formed  by  each  cell  class  are  presented  in  Supplementary 
Table  1 .  There  are  seven  cell  classes  modeled:  four  classes  of  exci¬ 
tatory  cells  including  layer  ll/lll  pyramidal  cells,  layer  IV  stellate 
cells,  layer  V  pyramidal  cells,  and  layer  VI  pyramidal  cells,  and 
three  classes  of  interneurons  including  basket  cells,  double  bouquet 
cells,  and  chandelier  cells.  Most  of  the  model  changes  described 
in  the  described  studies  involve  alterations  in  connection  num¬ 
bers  between  layer  ll/lll  pyramidal  cells,  one  of  the  known  robust 
horizontal  connections  systems  in  the  cortex  supporting  epileptic 
propagation  (Telfeian  and  Connors,  1998).  The  base  connection  for 
this  system,  6/2/3:2/3  =  178,  is  defined  as  the  number  of  layer  ll/lll 
pyramidal  cells  a  given  layer  ll/lll  pyramidal  cell  contacts  in  its 
axonal  distribution. 

The  model  in  general  illustrates  consistent  bursting  behavior, 
with  epochs  of  spontaneous  bursting  onset  and  cessation  given  a 
random  background  input  of  Poisson  based  charge  injection  to  1% 
of  the  cells  in  the  model.  This  is  an  effort  to  treat  the  underlying 
cortical  activity  as  input  from  neighboring  cortex,  with  the  model 
itself  treated  as  the  epileptic  focus  given  its  ability  to  produce  net¬ 
work  bursting  epochs.  The  synaptic  input  used  for  the  background 
was  not  periodic  in  nature.  Average  rates  for  these  Poisson  distri¬ 
butions  are  described  in  "Activity  changes  with  mean  background 
frequency”  section  and  Fig.  2.  The  synaptic  activations  used  for  the 


Statistics  and  analysis 

The  interbursting  phase  intervals  in  the  model  were  fit  with  a 
gamma  distribution  (Suffczynski  et  al.,  2005,  2006).  The  functional 
form  of  this  distribution  /(. . .)  is  given  by 


/(At)  =  (^r(a)r1  A  t""1  exp 


where  At  is  the  interbursting  phase  interval,  a  is  the  shape 
parameter,  p  is  the  scale  parameter,  and  r(...)  represents  the 
gamma  function.  Parameters  were  estimated  using  the  MATLAB 
function  gamfit  which  returns  maximum  likelihood  estimates  and 
95%  confidence  intervals  for  the  shape  and  scale  parameters.  A  non- 
parametric  test  of  randomness  was  used  to  attempt  to  establish 
whether  intervals  were  being  drawn  at  random  from  their  underly¬ 
ing  joint  distribution.  This  was  based  on  the  circular  definition  of 
the  lag-1  serial  correlation  coefficient  (Wald  and  Wolfowitz,  1943). 
p-Values  were  computed  under  the  assumption  of  asymptotic  nor¬ 
mality  of  the  test  statistic. 

A  method  for  the  detection  of  unstable  periodic  orbits  (of  period- 
1)  in  successive  interbursting  phase  intervals  was  applied  (So  et  al., 
1996,  1997).  A  period-1  orbit  is  a  fixed  point  of  the  nonlinear  map 
expressing  the  evolution  of  the  state  of  a  system,  iterated  a  single 
time  (Guckenheimer  and  Holmes,  1983).  Intervals  were  embedded 
in  a  two  dimensional  state  space  and  104  sets  of  transformed  inter¬ 
vals  were  obtained  after  randomization.  Dimensional  reduction  was 
instituted  using  circles  of  radius  9.4  ms  centered  along  the  diago¬ 
nal  of  the  state  space  (Le  Van  Quyen  et  al.,  1997).  One  hundred 
surrogates  were  produced  to  test  the  significance  of  peaks  which 
appeared  along  this  diagonal.  The  surrogates  were  generated  using 
the  amplitude  adjusted  Fourier  transform  algorithm  (Theiler  et  al., 
1992).  This  shuffles  the  original  sequence  of  interbursting  phase 
intervals,  maintaining  the  original  amplitude  distribution  of  the 
data  while  approximately  matching  its  Fourier  power  spectrum. 


Results 

A  total  of  1600  s  of  discontinuous  20-  and  30-s  data  seg¬ 
ments  were  obtained,  holding  the  base  connectivity  of  the 
layer  ll/lll  pyramidal  cells  to  layer  ll/lll  pyramidal  cells  at 
N2/3:2/3  =  178.  Only  the  random  number  seed  supplied  to  the 
background  input  generator  was  varied  for  each  of  these 
runs.  Additionally,  five  continuous  segments  of  data  were 
obtained  with  the  base  connectivity  set  at  N2/3:2/3  =  172  of 
lengths  320  s,  250  s  and  Nui-.ui  =  1 78  of  lengths  195  s,  1 40  s, 
and  208  s.  These  data  were  used  for  the  dynamical  sys¬ 
tems  analysis  presented  below.  In  addition  to  these  data, 
sixteen  20-s  runs  were  obtained  with  the  model  while  vary¬ 
ing  the  mean  background  input  frequency  at  the  base  level 
of  connectivity.  Five  20-s  runs  were  obtained  at  the  base 
connectivity  while  varying  the  temporal  pattern  of  the  back¬ 
ground  input,  and  ten  20-s  runs  were  obtained  with  a  fixed 
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«=  1.44  sec  I  =1.46  sec  1=1.48  sec 


Figure  1  (A)  Representative  connectivity  of  the  excitatory  cellular  component  in  a  given  modeled  minicolumn,  wiring  after 

(Douglas  and  Martin,  2004).  (B)  Three  dimensional  arrangement  of  the  16x16  array  of  minicolumns  in  space.  (C)  Representative 
snapshots  of  evolving  activity  over  0.02s  in  the  layer  ll/lll  pyramidal  cell  component.  Each  pixel  represents  one  cell,  color  coded 
proportionally  to  the  number  of  action  potentials  fired  in  bins  of  1/100  of  a  second. 
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Figure  2  Network  activity  produced  by  sequential  increases  in  the  mean  frequency  of  the  applied  background  activity  (background 
synaptic  input  provided  to  a  fixed  1%  set  of  the  modeled  cells,  summed  layer  ll/lll  pyramidal  cell  action  potentials  in  10  ms  bins). 
The  model  exhibits  a  transition  from  episodic  bursting  to  a  very  regular  bursting  behavior  driven  by  the  background  input. 


sequence  of  background  input  while  varying  the  overall  layer 
ll/lll  to  layer  ll/lll  connectivity  (N2/ 3:2/3). 


Activity  changes  with  mean  background  frequency 

These  experiments  were  performed  with  the  base  connec¬ 
tivity  of  the  layer  ll/lll  pyramidal  cell  system  set  at  the 
base  value  of  N2/ 3:2/3  =  178.  If  the  mean  frequency  of  the 
applied  background  synaptic  input  is  varied  from  0.25/s 
up  to  10/s,  several  patterns  of  activity  become  apparent 
(Fig.  2).  At  0.25/s,  only  the  low  level  set  of  activity  pro¬ 
duced  with  action  potential  production  by  the  background 
input  is  observed,  at  the  cells  where  the  input  takes  place. 
This  plot  essentially  demonstrates  the  Poisson-based  random 
network  activity  between  the  bursting  phases.  At  0.5/s,  spo¬ 
radic  bursting  activity  transmitted  to  the  network  as  a  whole 
can  be  observed,  with  long  quiescent  epochs.  As  the  applied 
mean  background  frequency  is  further  increased,  longer  and 
longer  periods  of  constant  bursting  activity  can  be  observed 
up  to  5/s.  After  this,  a  second  activity  transition  is  observed, 
in  which  the  activity  changes  from  continuous  bursting  into 
short  periods  of  very  large  amplitude  bursting  (in  terms  of 
numbers  of  active  neurons)  punctuated  by  brief  periods  of 
quiescence.  At  10/s  it  appears  to  dominate  the  activity.  This 
implies  a  saturation  mechanism  in  this  class  of  connected 
network,  which  comes  into  play  after  a  critical  percentage 
of  cells  are  excited  per  time  step.  This  saturation  behav¬ 
ior  is  again  an  intrinsic  property  of  the  fixed  network  being 
probed.  Additionally,  these  studies  imply  that  seizure  onset 
can  be  driven  by  neighboring  cortical  activity,  albeit  regu¬ 
lar  patterns  that  might  not  be  typical  of  random  background 
input  activity  utilized  here. 


Network  activity  altered  with  input  pattern 
changes 

Within  the  context  of  this  model,  it  is  possible  to  change  the 
random  pattern  of  connectivity  between  represented  cells, 
and  still  keep  the  total  number  of  connections  between  the 
various  cell  classes  constant.  By  varying  the  random  number 
seed  supplied  to  the  generator  distributing  the  connections, 
different  patterns  of  activity  can  be  demonstrated,  even 
with  the  same  application  of  underlying  cortical  activity 
applied  to  the  same  cells.  Examples  of  the  changes  in  activ¬ 
ity  are  presented  in  Fig.  3A.  The  pattern  produced  ranges 
from  almost  constant  bursting  throughout  the  20  s  exam¬ 
ined,  to  brief  periods  of  on  and  off  bursting.  Similarly,  the 
connectional  pattern  can  be  held  constant  along  with  the 
cells  in  which  the  background  activity  is  applied,  while  vary¬ 
ing  the  random  number  seed  responsible  for  producing  the 
order  in  which  background  pulses  are  injected  into  the  cells. 
This  produces  similar  alterations  in  network  activity  demon¬ 
strated  in  Fig.  3B,  and  can  include  several  time  scales  of 
bursting  epochs.  These  studies  imply  a  rich  dynamics  of 
stochastic  behavior  in  randomly  connected  neural  networks 
receiving  temporally  uncorrelated  background  input,  and 
again  point  toward  difficulties  in  predicting  when  the  burst¬ 
ing  phases  might  begin. 

Network  activity  is  very  sensitive  to  numbers  of 
excitatory  connections 

Finally,  changing  the  numbers  of  connections  in  this  net¬ 
work  model  can  produce  substantial  alterations  in  network 
behavior.  Fig.  4  presents  a  sequence  of  plots  of  the  layer 
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Figure  3  (A)  Network  activity  induced  by  varying  the  random  connectivity  pattern  between  cell  classes  (different  connectivity 

seeds,  N“£dn ,  for  the  random  number  generator).  Numbers  of  action  potentials  in  layer  ll/lll  pyramidal  cell  component,  10ms  time 
bins.  (B)  Network  activity  induced  by  varying  the  random  time  sequence  of  background  synaptic  input,  j.  In  these  experiments, 
all  cellular  connections  remain  fixed,  and  the  identity  of  the  cells  undergoing  background  synaptic  input  remain  fixed. 


ll/lll  pyramidal  cell  activity  (time-binned  action  potential 
numbers)  for  various  degrees  of  extracolumnar  connectivity 
between  the  layer  ll/lll  pyramidal  cells.  The  base  activ¬ 
ity  explored  in  this  manuscript  is  shown  in  the  plot  with 
total  number  of  connections  held  at  N2/3-.2n  =  178  (between 
extra-minicolumnar  layer  ll/lll  pyramidal  cells).  A  rapid 
reduction  in  network  bursting  is  shown  for  a  connection  num¬ 
ber  reduced  below  this,  and  almost  continuous  activity  is 
shown  for  connection  numbers  above  this.  Fig.  3A  and  B  data 
were  obtained  with  the  connectivity  set  at  N2/3:2/3  =  178. 
Only  the  random  pattern  of  connectivity  is  varied  in  Fig.  3A 
and  the  time  sequence  of  background  input  in  Fig.  3B.  These 
studies  were  performed  in  the  context  of  a  constant  aver¬ 
age  level  of  surrounding  background  input,  and  imply  the 
importance  of  internal  connectivity  in  the  development  of 
uncontrolled  bursting  of  the  network. 


Statistical  analysis  of  interburst  phase  intervals 

The  statistical  properties  of  the  interburst  phase  intervals 
for  five  continuous  runs  of  the  model  at  connectivities  of 
^2/3:2/3  =  178  and  172  were  analyzed.  This  was  motivated 
from  a  dynamical  systems  perspective,  where  periodicity  in 
sequential  intervals  was  sought  for.  Fig.  5A  displays  a  his¬ 
togram  approximation  to  the  probability  density  function 
for  interbursting  phase  intervals  for  Continuous  Run  1 ,  which 
consisted  of  163  intervals  collected  from  a  320  s  run  of  the 
simulation.  Fig.  5D  displays  the  same  histogram  approxi¬ 
mation  for  Continuous  Run  4,  with  a  total  of  58  intervals 
collected  from  a  140s  run  of  the  simulation.  Gamma  distri¬ 
butions  were  used  to  fit  these  densities  (Fig.  5A  and  D  (blue 
traces))  (see  "Materials  and  methods"  section;  Suffczynski 
et  al.,  2005,  2006).  In  the  case  of  Continuous  Run  1 
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Figure  4  Changing  the  absolute  connectivity  in  the  layer  ll/lll  pyramidal  cell  component  (number  of  layer  ll/lll  pyramidal  cells 
contacted  by  a  given  layer  ll/lll  pyramidal  cell,  N2/3-.2b)  in  the  model  produces  alterations  in  the  network  bursting  behavior.  At 
very  low  absolute  connectivity  {Ni/i-.in  =  110)  network  bursting  is  brief  and  isolated,  while  at  higher  levels  of  absolute  connectivity 
periods  of  constant  bursting  can  be  observed. 


(Fig.  5A)  we  found  a  =  6.09  (4.61—8.04)  for  the  shape  param¬ 
eter  of  the  distribution,  and  ,6  =  0.26  (0.20—0.35)  for  the 
scale  parameter  (95%  confidence  intervals  in  parentheses). 
For  Continuous  Run  4  (Fig.  5D),  the  corresponding  values 
were  or  =  7. 1 4  (5.00—10.19),  and  ^  =  0.23  (0.16—0.33)  (simi¬ 
lar  results  were  obtained  for  3  further  continuous  runs  and 
one  discontinuous  run  —  results  not  shown).  In  accord  with 
the  interpretation  of  Suffczynski  et  al.  (2005,  2006),  the  fact 
that  the  shape  parameter  a  was  larger  than  one  in  all  runs, 
suggests  the  potential  presence  of  periodicity  in  the  genera¬ 
tion  of  bursting  epochs.  To  probe  this  link  further,  additional 
statistical  tests  were  performed  as  described  below. 

To  ascertain  how  intervals  were  being  drawn  from  their 
underlying  joint  distribution,  we  applied  a  non-parametric 
test  of  randomness  to  the  interbursting  phase  intervals  (Wald 
and  Wolfowitz,  1943).  In  the  case  of  Continuous  Run  1,  we 
found  that  one  cannot  reject  the  null  hypothesis  of  ran¬ 
domness  at  the  5%  significance  level  (p-value  0.93).  For 
Continuous  Run  4,  we  found  the  same  conclusion  at  the  5% 
significance  level  (p-value  0.29).  This  conclusion  was  also 
borne  out  for  the  remaining  three  continuous  runs. 

A  method  for  the  detection  of  unstable  periodic  orbits 
(of  period-1)  was  then  applied  to  test  for  the  presence 
of  deterministic  dynamics  in  the  generation  of  interburst¬ 
ing  phase  intervals  (So  et  al.,  1996,  1997).  This  method 
institutes  a  transformation  of  the  sequence  of  intervals  such 


that  the  transformed  sequence  is  clustered  around  locations 
of  potential  periodic  orbits.  One  can  compare  the  peaks  of 
these  clusters  to  those  generated  by  surrogate  data  (Theiler 
et  al.,  1992),  to  compute  the  statistical  significance  of  the 
peaks,  and  thereby  ascertain  the  potential  existence  of  peri¬ 
odic  orbits  in  the  data.  Fig.  5C  and  F  shows  the  peaks  of  the 
clusters  were  not  significantly  greater  than  those  generated 
by  surrogate  data  (see  caption),  and  so  no  period-1  orbits 
were  detected  for  either  run  (nor  for  the  remaining  three 
continuous  runs),  at  the  limit  of  detection  in  the  current 
data  set. 


Discussion 

Our  results  demonstrate  that  while  holding  the  mean  prop¬ 
erties  of  the  network  stable  (mean  connectivity  numbers, 
mean  background  excitation  rates),  very  rich  and  strikingly 
different  dynamics  are  produced  by  changing  the  model 
details.  Epileptogenic  behavior  can  be  created  in  these  net¬ 
works,  as  described  above,  by  changes  in  the  random  pattern 
of  connectivity,  while  holding  fixed  the  intrinsic  active  or 
passive  membrane  properties  in  the  constituent  neurons. 
Such  changes  in  connectivity  could  be  analogous  to  changes 
in  underlying  connectivity  that  might  occur  following 
cerebral  insults,  or  repetitive  seizures.  Similar  modeling 
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Figure  5  Statistical  analysis  of  the  interbursting  phase  intervals.  (A  and  D)  Interval  histograms  are  fit  with  gamma  distributions 
revealing  shape  parameters  a>  1  in  the  case  of  Continuous  Run  1  (A)  and  Continuous  Run  4  (D).  (B  and  E)  Two  dimensional  delay 
embedding  of  sequential  intervals  (Spearman  correlation  coefficients  are  not  significantly  different  from  zero)  for  Continuous  Runs 
1  and  4,  respectively.  (C  and  F)  Testing  the  significance  of  potential  period-1  orbits  detected  through  a  dynamical  analysis  of  the 
sequence  of  points  displayed  in  (B)  and  (E),  respectively.  The  y-axis  represents  the  fraction  of  surrogate  (shuffled,  see  "Materials 
and  methods”  section)  sequences  with  a  maximal  deviation  from  the  mean  surrogate  result  of  greater  than  W  (So  et  al.,  1996, 
1997)  (104  random  matrices  and  102  surrogates  were  used).  The  horizontal  dotted  (red)  line  displays  the  maximal  deviation  for  the 
simulation  data.  Since  there  exists  a  significant  fraction  of  surrogates  with  deviation  greater  than  that  for  the  simulation  data  (for 
both  C  -  20%  and  F  -  10%),  neither  plot  displays  convincing  evidence  of  the  existence  of  a  period-1  orbit.  (For  interpretation  of 
the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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work  has  been  performed  on  changes  in  connectivity  with 
resultant  epileptic  effects  in  hippocampus  (Morgan  and 
Soltesz,  2008;  Dyhrfjeld-Johnsenetal.,  2007).  In  our  results, 
local  increases  or  decreases  in  connectivity  in  the  model  may 
alter  the  network  in  a  similar  manner  (a  tempting  compar¬ 
ison  which  must  be  tempered  by  our  lack  of  understanding 
of  the  functional  significance  of  new  or  absent  connections, 
see  Sallin  et  al. ,  1995;  Dinocourt  et  al. ,  2003;  Marco  et  al., 
1997;  Dudek  and  Sutula,  2007;  Magloscky,  2010).  In  the 
model  above,  the  external  background  activity  can  be  held 
constant  or  changed  in  a  variety  of  time-frequency  manip¬ 
ulations.  Various  simulations  can  be  created  where  seizures 
occur  rarely  or  very  frequently. 

Similarly,  human  epileptic  seizures  are  episodic,  tran¬ 
sient  events.  Whether  epileptic  seizures  are  random  events 
is  not  clear,  but  times  of  ictal  onset  can  behave  as  a  random 
process  (Suffczynski  et  al.,  2006).  In  some  patients  there 
are  no  identifiable  contributing  factors,  in  other  patients 
such  conditions  such  as  sleep  deprivation  may  increase  the 
chance  of  seizure  occurrence  and  in  still  other  patients 
seizures  can  be  provoked  by  specific  stimuli  (e.g.  hyperven¬ 
tilation,  intermittent  photic  stimulation;  Lu  et  al.,  2008; 
Vinogradova  et  al.,  2009;  Kaplan  et  al.,  2009).  In  these 
examples  there  may  be  resulting  alterations  of  background 
activity  within  and  outside  the  epileptogenic  zone  (not  evi¬ 
dent  from  the  EEG)  that  make  it  more  likely  that  a  seizure 
may  arise  from  the  existing  focal  epileptogenic  network. 
While  these  influences  might  create  a  state  where  seizures 
are  more  likely  to  occur,  this  facilitatory  state  should  be  dis¬ 
tinguished  from  the  presence  (or  absence)  of  a  preictal  state 
in  unprovoked  seizures. 

Random  behavior  of  network 

One  interesting  aspect  of  this  model  is  its  ability  to  demon¬ 
strate  both  spontaneous  periods  of  bursting  activity  as 
well  as  self-termination  of  the  bursting.  As  illustrated  in 
Figs.  2—5,  these  periods  of  bursting  can  be  quite  variable  in 
their  length.  Throughout  the  bursting  and  quiescent  phases, 
the  distributed  background  activity  is  constantly  active, 
affecting  approximately  1%  of  the  total  cells.  We  believe 
this  model  may  represent  the  type  of  network  behavior 
described  by  Lopes  da  Silva  et  al.  (2003),  in  which  epileptic 
activity  within  the  network  cannot  be  predicted  from  the 
interictal  state.  The  non-parametric  and  nonlinear  dynami¬ 
cal  analyses  described  in  the  "Results”  section  support  this, 
however  inference  from  this  is  somewhat  limited  given  the 
finite  size  of  the  data  set,  the  small  region  of  modeled  area, 
and  the  gamma  distribution  fits  described  earlier. 

This  obviously  may  not  be  true  of  all  types  of  clini¬ 
cal  and  model  epileptiform  behavior.  For  example  Osorio 
et  al.  (2009,  2010)  show  that  pharmacoresistant  seizures 
tend  to  cluster,  and  may  have  an  inherent  self-triggering 
capacity.  This  might  make  a  prediction  algorithm  possible 
to  implement  in  a  useful  fashion.  Others  (Suffczynski  et  al., 
2005,  2006),  however  demonstrate  that  seizure  onset  can 
be  described  in  both  experimental  and  model  data  as  a  ran¬ 
dom  walk  process,  with  possibly  a  deterministic  mechanism 
ascribed  to  seizure  termination.  The  random  onset  nature 
would  (in  a  bistable  network  with  Poisson  transitions)  be 
difficult  to  predict. 


In  the  case  of  our  model  the  underlying  fluctuations  lead¬ 
ing  to  seizure  onset  are  the  random  background  activity  we 
have  imposed  on  the  network.  The  properties  of  the  network 
connectivity  then  support  the  bursting  frequency  observed 
(Anderson  et  al.,  2007).  These  types  of  grossly  synchronized 
bursting  states  in  the  context  of  a  neural  network  have 
been  studied  extensively  by  Kowalski  et  al.  (1992).  They  are 
truly  pathological  in  the  sense  that  they  would  block  or  con¬ 
found  any  information  flow  through  this  network.  It  is  also  a 
network-generated  state,  and  can  be  stopped  by  eliminating 
synaptic  transmission  (Kowalski  et  al.,  1992;  Keefer  et  al., 
2001;  Rhoades  and  Gross,  1994). 

Many  limiting  cases  of  the  gamma  distribution  have  phys¬ 
ical  interpretations  that  might  make  it  easier  to  understand 
spiking  data  from  cortex  when  used  for  fitting  (Papoulis, 
1984;  Suffczynski  et  al.,  2005).  When  the  shape  parameter, 
a,  in  the  gamma  distribution  is  an  integer,  the  distribution 
is  known  as  an  Erlang  distribution  and  represents  the  prob¬ 
ability  distribution  of  the  waiting  time  until  the  ath  event 
from  a  sampled  Poisson  process  with  characteristic  time  p. 
This  might  be  comparable  to  the  distribution  of  the  number 
of  spiking  or  underlying  synaptic  events  required  to  trigger 
the  network  bursting  behavior.  One  could  envision  trying  to 
extract  the  integer  value  of  a  from  either  computational 
or  experimental  data.  Similarly,  the  Maxwell— Boltzmann 
distribution  can  be  related  to  the  gamma  distribution 
under  certain  restrictions  on  the  gamma  scale  and  shape 
parameters,  implying  the  possibility  of  extracting  almost 
thermal-like  or  statistical  mechanical  interpretations  of  the 
network  activity  (Hegyi,  1996). 

Limitations  of  the  model 

Our  plots  demonstrate  gross  summed  numbers  of  time- 
binned  action  potentials  in  the  model  for  given  neuron 
classes.  The  interictal  activity  is  a  random  Poisson  input  to 
1%  of  the  cells  in  the  model  and  is  demonstrated  in  the  plots, 
particularly  in  Fig.  2,  Panel  1,  inset.  We  chose  to  impose 
this  random  interictal  behavior  on  the  model  to  demon¬ 
strate  that  these  fluctuations  can  produce  very  coherent 
synchronous  oscillations  in  an  unpredictable  fashion.  Flow- 
ever,  the  surrounding  input  might  not  have  to  be  completely 
random  to  bring  about  the  same  effect.  Epilepsies  involving 
specific  stimuli  might  require  a  coherent  surrounding  input 
to  give  rise  to  the  seizure  (Lu  et  al.,  2008;  Vinogradova  et  al., 
2009;  Kaplan  et  al.,  2009).  A  more  realistic  technique  would 
be  to  treat  the  interictal  background  as  a  log-normal  process 
which  does  have  some  support  in  the  literature  (see  Farkhooi 
et  al.,  2009;  Waters  and  Helmchen,  2006).  Newer  recording 
methods  from  invasively  monitored  epilepsy  patients  might 
help  determine  what  patterns  of  background  activity  are 
causative  (Truccolo  et  al.,  2011;  van  Gompel  et  al.,  2008). 
Our  intent  was  to  represent  the  resting  interictal  cellular 
activity  as  fundamentally  sparse  with  a  random  component. 
This  was  most  easily  implemented  as  a  low  frequency  Poisson 
process. 

It  is  possible  to  view  the  single  synaptic  input  (driven 
by  the  background  source)  as  representing  several  weaker 
but  synchronized  inputs.  This  represents  a  limitation  to  this 
modeling  approach,  a  limitation  that  in  large  degree  could 
be  corrected  with  more  elaborate  multicompartment  and 
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synaptic  representations  of  the  cells  to  make  the  multiple 
weaker  inputs  more  independent.  Additionally  in  this  lim¬ 
ited  data  set,  we  are  unable  to  say  much  about  time  epochs 
larger  than  several  hundreds  of  seconds  (our  largest  con¬ 
tinuous  simulation  being  320s).  This  is  clearly  a  limitation 
in  this  technique  and  future  computational  work,  including 
efforts  in  our  laboratory  will  explore  longer  time  intervals 
of  ictal  and  interictal  behavior,  and  possibly  push  the  detec¬ 
tion  limit  for  predictable  activity  lower  (or  detect  it  more 
accurately).  Models  such  as  this  particular  rigid  crystalline 
arrangement  of  neurons  with  many  fixed  cellular  properties 
in  some  sense  have  less  inherent  "randomness”  than  real 
neocortex.  This  work  was  primarily  meant  to  spark  inter¬ 
est  in  a  possible  mechanism  for  the  difficulties  inherent  in 
seizure  prediction,  but  by  no  means  should  it  be  interpreted 
too  literally. 

Implications  for  seizure  prediction 

The  purpose  of  the  model  presented  here  is  not  to  judge  the 
effectiveness  of  seizure  prediction,  but  rather  to  present  a 
plausible,  alternative  hypothesis  for  partial  seizure  occur¬ 
rence  that  could  explain  situations  where  seizure  prediction 
may  not  be  possible.  Indeed,  it  is  conceptually  attractive  to 
consider  that,  just  as  partial  seizures  may  result  from  vari¬ 
ous  pathologies  and  mechanisms,  that  some  partial  seizures 
may  not  be  reliably  predicted.  It  is  beyond  the  scope  of 
this  discussion  to  address  the  various  methods  being  used 
in  attempts  to  predict  epileptic  seizures.  It  is  always  impor¬ 
tant  to  differentiate  true  seizure  prediction  from  improved 
seizure  detection.  Other  commentaries  and  reviews  address 
these  methods  and  include  discussions  of  the  challenges 
and  frustrations  to  date  in  routine  seizure  prediction  even 
with  intracranial  electrodes  (Estellar  et  al.,  2001;  Litt  and 
Echauz,  2002;  Sackellares  et  al.,  2006;  Haas  et  al.,  2007; 
Osorio  et  al.,  2001). 

This  study  focused  on  the  interval  to  the  time  of  the 
next  "seizure”  or  busting  phase  in  the  model.  Our  interest 
was  in  the  occurrences  of  the  transitions  from  the  quiescent 
or  background  state  into  the  pathologic  state,  since  that 
is  what  most  seizure  prediction  algorithms  are  optimized 
for.  There  is  fairly  strong  evidence,  certainly  in  the  case  of 
complex  partial  seizures  in  temporal  lobe  epilepsy,  that  the 
length  of  clinical  seizures  can  be  fairly  uniform  in  a  given 
patient  (see  for  instance  Afra  et  al.,  2008).  The  time  inter¬ 
val  durations  of  the  seizures  themselves  may  also  prove  to 
have  to  predictive  guidance  as  well  and  should  be  explored 
in  the  future  in  modeling  efforts.  This  may  be  more  useful 
in  the  case  of  neocortical  epilepsy  with  its  rapid  spread  and 
possible  involvement  of  larger  regions  of  tissue. 

Additionally,  this  model  can  incorporate  incremental 
changes  in  connectivity  in  the  epileptogenic  zone,  changes 
that  could  be  a  model  for  progressive  epileptogenesis  (e.g. 
sprouting).  This  type  of  model  also  provides  data  that  is 
comparable  to  clinical  data  from  epilepsy  patients.  The 
simulated  network  activity  is  taken  from  a  small  region  of 
modeled  cortex  comparable  in  size  to  the  surface  area  under 
a  typical  subdural  grid  electrode,  and  makes  comparisons 
between  modeling  efforts  and  clinical  data  easy  to  perform 
(Anderson  et  al.,  2007,  2009;  Kudela  et  al.,  1997,  2003a, b; 
Franaszczuk  et  al.,  2005).  Indeed  the  major  advantage  of 


neuronal  network  modeling  is  the  ability  to  simultaneously 
monitor  activity  in  all  of  the  network  neurons  under  given 
experimental  conditions,  something  not  possible  with  bio¬ 
logical  systems,  even  with  sophisticated  recording  arrays. 
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